# libs ----
library(sf)
library(geomander)
library(tidyverse)
library(redist)
library(ppmf)
library(grid)
library(patchwork)
library(wacolors)

# consts ----

state <- 'PA'
# get data ----

# data ----
ppmf12_path <- "../../data-raw/ppmf_12.csv"
ppmf4_path <- "../../data-raw/ppmf_04.csv"

ppmf12 <- read_ppmf(state, ppmf12_path)
ppmf4 <- read_ppmf(state, ppmf4_path)

ppmf12 <- ppmf12 %>% add_geoid() %>% agg() %>% breakdown_geoid()
colnames(ppmf12) <- paste("v12", colnames(ppmf12), sep = "_")
ppmf12 <- ppmf12 %>% rename(GEOID = v12_GEOID)

ppmf4 <- ppmf4 %>% add_geoid() %>% agg() %>% breakdown_geoid()
colnames(ppmf4) <- paste("v4", colnames(ppmf4), sep = "_")
ppmf4 <- ppmf4 %>% rename(GEOID = v4_GEOID)

# comparison ----
census <- create_block_table(state = state)

# all joined ----
block <- census %>%
  left_join(ppmf12, by = 'GEOID') %>%
  left_join(ppmf4, by = 'GEOID')

# and remove duplicates
block <- block %>% select(-contains('.'))

# and set missing block pop/vap to 0
block[is.na(block)] <- 0

# add block_group back (dropped by contains('.'))
block <- block %>% breakdown_geoid()
congress <- tigris::congressional_districts(state == "PA")
congress <- congress %>% filter(STATEFP == 42)

# match ----
block_congress_match <- geo_match(from = block,
                                to = congress,
                                method = 'centroid')

# Add to data ----
block <- block %>%
  mutate(cd = congress$CD116FP[block_congress_match])

# clean data ----
dist <- block %>%
  st_drop_geometry() %>%
  group_by(cd) %>%
  summarize(pop = sum(pop),
            pop4 = sum(v4_pop),
            pop12 = sum(v12_pop))

# parity ----
(en_p  <- redist.parity(as.numeric(dist$cd), dist$pop))
(v4_p  <- redist.parity(as.numeric(dist$cd), dist$pop4))
(v12_p <- redist.parity(as.numeric(dist$cd), dist$pop12))

# pop diff ----
en_pop  <- max(abs(dist$pop - mean(dist$pop)))
v4_pop  <- max(abs(dist$pop4 - mean(dist$pop4)))
v12_pop <- max(abs(dist$pop12 - mean(dist$pop12)))

pa_parity <- tibble(names = c("Census", "DAS-4.5", "DAS-12.2"),
       parity = c(en_p, v4_p, v12_p),
       population = c(en_pop, v4_pop, v12_pop))

# for use in paper ----
# saveRDS(pa_parity, "../../data/PA/pa_parity.rds")
